import matplotlib.pyplot as plt     #画图包
import numpy as np                  #科学计算包
from scipy import integrate, linalg, signal  #数值积分计算
from scipy.sparse.linalg import eigs #计算特征向量
from scipy.integrate import odeint  #计算微分方程

def LinerTransform(data):
    c=0  #线性偏差常数，用于主函数预测
    x0=np.array(data,float) #将数据转化为浮点数，提高计算精度
    len=x0.shape[0]    #获取数据长度
    #计算级比的上下限
    low=np.exp(-2/(len+1))  
    high=np.exp(2/(len+2))
    #获取数据中最大值、最小值
    x_max=x0.max()
    x_min=x0.min()

    #计算级比,同时获取级比中的最大值和最小值
    JiBi=x0[0:-1]/x0[1:]
    JiBi_max=JiBi.max()
    JiBi_min=JiBi.min()

    #将级比迭代，直到满足达到要求，在范围之内
    while(JiBi_max>high or JiBi_min<low):
        x0=x0+x_min  #最小值作为偏差，将数据加上一个偏差
        c=c+x_min    #更新线性偏差
        #更新数据中的最大值最小值，便于后续迭代
        x_max=x0.max()  
        x_min=x0.min()
        #更新级比
        JiBi=x0[0:-1]/x0[1:]
        JiBi_max=JiBi.max()
        JiBi_min=JiBi.min()

    return c  #返回线性偏差，与原数据相加符合灰色预测要求

#主函数/灰色预测模型
def GreyModel(data,predition):
    x0=np.array(data,float) #将数据转化为浮点数，提高计算精度(同上线性偏差函数)
    len=x0.shape[0]    #获取数据长度
    x1=np.cumsum(x0)   #累加生成序列、生成(1-AGO)序列
    B = np.array([-0.5*(x1[0:-1]+x1[1:]), np.ones(len-1)]).T #构建B矩阵
    Y=x0[1:]   #构建Y向量，便于后续计算最小二乘解析法中的参数  
    #最小二乘估计参数
    a_hat=linalg.lstsq(B,Y)[0] 

    #定义差分方程，用于计算微分方程的数值解
    def Diff(y,t,a,b):
        return np.array(-a*y+b)
    
    t=np.arange(len+predition)  #定义时间向量
    solve=odeint(Diff,x0[0],t,args=(a_hat[0],a_hat[1])) #利用odeint函数求解微分方程的数值解
    solve=solve.squeeze()  #压缩其中的零向量解
    result = np.hstack((x0[0], np.diff(solve)))  #计算差分结果，即就是我们需要的预测值
    return result

#原始数据,即为2014~2023每年杭州的平均房价
data=np.array([17730,16303, 18892, 26365, 28856, 27944, 26866, 30742, 34356,32549])
#线性偏差c
c=LinerTransform(data)
#提高主函数计算出来的预测值,因为加上了线性偏差，所以计算出来的结果需要减去线性偏差
data_hat=GreyModel(data+c, 2)-c

#画图
t1 = range(2014, 2024)
t2 = range(2014, 2026)
plt.plot(t1, data, color='r', linestyle="-", marker='*', label='True')
plt.plot(t2, data_hat, color='b', linestyle="--", marker='.', label="Predict")
plt.legend(loc='upper right')
plt.xlabel('years')
plt.ylabel('yuan/square meter')
plt.title('Prediction by GM(1,1)')
plt.gcf().set_facecolor('lightgrey')
plt.show()

#以下代码为两个检验函数代码
# 计算标准差
def StandardDeviation(x, y):
    x=np.array([x, np.ones(x.size)]).T
    u=linalg.lstsq(x, y)    #最小二乘法拟合线性模型
    residuals=y-np.dot(x, u[0]) 
    std_dev=np.std(residuals) #计算残差的标准差
    return std_dev

#相对残差
def RelativeResidual(data,pre_data):
    return np.abs((data-pre_data)/data)

#调用函数计算
relative_residual=RelativeResidual(data,data_hat[:len(data)])
standard=StandardDeviation(data_hat[0:data.size], data)
print("标准差为：",standard)
print("相对残差：",relative_residual)

